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A new registration algorithm based on 
Newton-Raphson iteration is proposed to 
align images with rigid body transforma- 
tion. A set of transformation parameters 
consisting of translation in x and v and 
rotation angle around z is calculated by 
optimizing a specified similarity metric 
using the Newton-Raphson method. This 
algorithm has been tested by registering 
and correlating pairs of topography 
measurements of nominally identical NIST 
Standard Reference Material (SRM 2461) 
standard cartridge cases, and very good 
registration accuracy has been obtained. 
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1, Introduction 

Image registration is a process of determining the 
point-by-point correspondence between two or more 
images taken of a scene at different times, by different 
sensors, or from different points of view. The images 
are aligned with one another so that differences 
between them can be detected. The method has applica- 
tions in many areas such as remote sensing, medical 
imaging, computer vision, image similarity measure- 
ment, etc. Based on different criteria, registration 
methods and applications can be classified differently. 
Over the years, a broad range of techniques has been 
developed for registration of various types of data and 
images [1,2]. 

The parameters that make up the registration trans- 
formation can either be computed directly, if possible, 
or be determined by finding an optimum of some 
function defined on the parameter space. In the latter 
case, the transformation parameters can be incorporat- 
ed into a standard mathematical function. This fimction 
attempts to quantify the similarity between two images 
given a certain transformation [2]. The goal of image 
registration is to find the optimum parameters to make 



the function reach its extremum. Various optimization 
techniques can be used to accomplish this. As a well- 
known numerical analysis method for finding succes- 
sively better approximations to the zeroes of a real- 
valued function, the Newton-Raphson method is one of 
the widely used optimization techniques [3]. With this 
method applied, Lucas and Kanade [4] presented a 
registration technique making use of the spatial intensi- 
ty gradient of images, Vendroux and Knauus [5] 
optimized a digital image correlation (DIG) coefficient 
to calculate surface displacements and displacement 
gradients of a sample under deformation, and Gole- 
Rhodes et al. [6] calculated a mutual information 
metric, which may be especially useful for registration 
of multi-modality images (i.e., images taken by differ- 
ent imaging methods) [7]. 

In response to a request from the Bureau of Alcohol, 
Tobacco, Firearms, and Explosives (ATF) to develop 
physical standards for bullets and casings for use in the 
validation of image analysis systems for ballistics 
identification, the National Institute of Standards and 
Technology (NIST) has developed the Standard 
Reference Material (SRM 2460) standard bullets and 
2461 standard cartridge cases [8]. The similarity 



Volume 115, Number 1, January-February 2010 

Journal of Research of the National Institute of Standards and Technology 



between the master standard bullet or cartridge case 
and their replicas, or between different replicas needs to 
be determined. This calculation requires an explicit 
similarity metric. A parameter called the signature 
difference £>, has been proposed by NIST to quantify 
the similarity for both 2D and 3D ballistics signatures 
[9]. Accordingly, an appropriate image registration 
method is needed to align two casing images before the 
similarity between them is calculated. For the practical 
requirement of determining the reproducibility of the 
topography signatures of bullets or cases, a new image 
registration method based on Newton-Raphson itera- 
tion is developed, and this is directly associated with 
the signature difference parameter In this paper, we 
introduce the mathematical algorithm in Sec. 2, we 
define measures of similarity in Sec. 3, and we describe 
the experimental tests in Sec. 4. 



and the transformation array needs to be solved, a 
registration technique is required. The object of 
registration is to determine the rigid body transforma- 
tion P that aligns the two acquired images, / and g, as 
well as possible. 

Let Sp represent any single point in the reference 
image / with the lateral coordinates (x, y) and P an 
estimation value of the real transformation P. For each 
estimation P, Sp has a corresponding point in the com- 
parison image, which is denoted as Sp. The intens- 
ity values of Sp and Sp can be denoted as f{Sp) and 
g{Sp, P), or equivalently as/(x, y) and g(x , y ) . We use 
a parameter C to quantify how well the two images 
match. It is defined by the ratio of sum of squared 
intensity differences between the two images to the 
sum of squared intensities of the original image, which 
is taken as the reference. See Eq. (3). 



2, The Mathematical Algorithm 

In image processing and computer vision, rigid body 
transformation can be decomposed into a rotation and a 
translation. We assume that the topographic images are 
embedded in a common coordinate system. If the trans- 
formation is in a 2D plane, which is a common simpli- 
fication in many optical image processing and measure- 
ment situations, the coordinate position (x, y) of any 
point inside the object being transformed will be 
relocated to (x, y ) according to 
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We assume that/and g, which are called respective- 
ly the reference image and the comparison image in this 
paper, are images acquired before and after the object is 
physically transformed and that the transformation is 
that of a rigid body. Let P represent the array of all the 
parameters of the transformation. In the case of the 2D 
rigid body transformation as defined by Eq. (2), P is 
composed of three transformation parameters, T^, Ty, 
and 9. For each transformation P, image g will be 
uniquely determined, and vice versa. When the refer- 
ence image, /, and the comparison image, g, are given 
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C is a function of P, the array of estimated transfor- 
mation parameters. When the correct P that transforms 
the object is substituted into Eq. (3), C will reach its 
minimum. The minimum should be zero for the ideal 
condition that noise and error are not considered. 
Conversely, an estimation P that minimizes C is the 
correct one. 

The pixels of both acquired reference and compari- 
son images are at grid positions. However, the coordi- 
nate positions (x, y ) of points in the comparison image 
mapped by P may not be at grid points. Therefore, in 
order to calculate the value of C, an interpolation 
process is required to define the function value g at 
non-grid positions. In order to obtain good piecewise 
continuity and smoothness of the interpolation function 
and avoid oscillation of the function, a bi-cubic spline 
interpolation [10] is applied using common methods 
[11]' for each choice of the array P. 

Now the goal to find an estimation array P to mini- 
mize the function C is an optimization problem. To find 
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the minimum of C, the gradient of C with respect to 
variation in the parameters of P must be zero. If T^, Ty, 
and Gave denoted aspi(i = 1,2,3), we will get 



VC = 



-2 
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where 
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Based on Eq. (4), it can be expanded as 
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The nonlinear Eq. (4) can be solved using the 
Newton-Raphson iteration method [3]. Then, we obtain 



p*H-i = p* _ VC(P*)-[VVC(P*)]" 



(6) 



where P * is the A:* iteration value and P *^' is 
the approximation value after the (k+ I) iteration. 
VVC(P) is the Hessian matrix [12]. 



If 
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are denoted as if^, they will be 
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Eq. (9) continued 
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When P is close to the exact solution, 
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(10) 



so that the first sum item of Eq. (8) can be ignored and 
we will get 
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3, Definitions and Measures of 
Registration 

In order to quantify the reproducibility of the topog- 
raphy of the SRM bullets and cartridge cases, we use 
two parameters that compare pairs of 3D images or 2D 
profiles, the cross-correlation function maximum 
{CCF.^ and a new parameter, the signature difference 
D, [9]. Given two casing signature images, denoted as/ 
and g following the registration algorithm description 
in Sec. 2, we take / as reference and transform the 
position of g using three parameters, horizontal trans- 
lation T^, vertical translation T^, and rotation angle 9. 
For a certain group of variants {T^, T^, 9), a correspon- 
ding image g' transformed from g, and the overlapping 
subarea of/and g', which are denoted as / and g'„ are 



determined. In the image processing application, 
the normalized cross-correlation function will be 
calculated as 



CCF{f,g,T,,T^„9)=CCF{f^,,g:)= 




(12) 



where < • , • > is the inner product and | • | is the 

i^-norm. 

The maximum value CCF,„^ occurs when the two 
correlated surface topographies are registered at their 
best matched position. Although the CCF,„^ can be 
used for casing signature comparison, it is not sensitive 
to z-scale differences. Assuming two compared signa- 
tures have the same shape but different amplitude 
scales, their CCF„,^. is still 100%. Therefore another 
parameter called the signature difference, D„ which is 
highly correlated with the CCF.^^^ but directly quanti- 
fies bullet profile signature difference, has also been 
proposed [9]. 

The parameter Z), is calculated in the following way. 
After the registration algorithm is performed, a trans- 
formed image g' which is at the best position to match 
the reference image/is determined. At this position, the 
overlapping areas in image /and image g' are denoted 
as/o and g'^o, respectively. The signature difference D, 
is defined as a ratio of the area-based mean-square 
roughness Sq^ of the topography difference g'j,o-/,o 
and the mean-square roughness of the reference topog- 
raphy / . The parameter Z), therefore casts the mean 
square difference as a fraction of the mean square 
roughness of the reference topography, so for good reg- 
istration of similar surfaces the value of £>, should be 
close to zero 



D, 



Sq'{gl,o-f,,o) 

Sq'iAo) 



(13) 



For a set v having n elements, the definition of the 
mean-square roughness Sq^ is 



Sq\v)-. 






(14) 



Equation (13) and Eq. (3) have the same physical 
meaning and they can be shown to be equivalent to one 
another More specifically, the parameter Z), is calculat- 
ed after the transform parameter P is solved and 
two images are aligned, whereas P is solved through 
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minimizing C, a parameter directly equivalent to D, . As 
a result, our registration algorithm, which searches a set 
of optimum transformation parameters to minimize the 
function C in Eq. (3), ensures the consistency between 
the image registration and the calculation of D,. It 
avoids possible error caused by inconsistent criteria 
between registration and similarity measurement. 



4, Experimental Data and Results 

The algorithm has been tested by registering three 
pairs of topographic images of standard cartridge cases 
to determine how similar they are. In all three cases 
very good correlation is obtained. An example is given 
as follows. Images shown in Fig. 1(a) and 1(b) are 
portions of topographic images of breech face impres- 
sions acquired from two different SRM cartridge cases 
by a commercial confocal microscope. They are filtered 
using a Gaussian filter [13] with 0.08 mm long-wave- 
length cutoff to remove the low frequency components 
of the original casing surface, i.e., the curvature and 
waviness. The filter process accentuates the individual 
characteristics impressed by the firearm. This system 
originally used an open source code to register different 
breech face or firing pin images before the similarity 
calculation [14]. However, during our recent test for 
breech face signature measurements, the software 
sometimes failed to register two correlated images at 
their best matched position. In this example, the CCF„,^^ 
calculated between two images after registration gave 
an extremely low value, CCF,„^j= 5.5 %. With the 
image in Fig. 1(a) selected as the reference, the image 
in Fig. 1(b) does not show any positional transform 



to match the reference after the registration process, 
indicating that the open source registration process did 
not work for this case. This error does not happen when 
applying the new registration algorithm described in 
Sec. 2. By applying this algorithm, an optimal array of 
transformation parameters P is determined. Then the 
comparison image is transformed back to match the ref- 
erence image using the array P as shown in Fig. 1(c). 
The resulting similarity metrics are CCF^^^ = 97.9 % 
and £>, = 4.2 %. Taking the center point of the image 
as coordinate origin, the calculated transformation 
parameters are T,. = 16.3 pixels, 7^=32.9 pixels and 
9=2.9°. 

All NIST standard cartridge cases are electroformed 
replicas from the same master casing provided by the 
ATF. High agreement and reproducibility between the 
casing signatures of different standard casings have 
been verified previously [15]. In this experiment, a 
0.08 mm long-wavelength cutoff Gaussiein filter has 
been adopted. This filter procedure emphasizes shorter 
surface spatial wavelengths than the 0.25 mm cutoff 
mainly used in our previous experiments and eliminates 
more components considered as low frequency wavi- 
ness of the topography. With this procedure, the exper- 
imental results still show that the proposed registration 
algorithm is an effective method for establishing the 
similarity of NIST standard cartridge cases and is capa- 
ble of registering surfaces with sub-pixel precision. 



5. Discussion 

The automatic algorithm proposed in this paper 
solves the 2D image registration problem for rigid body 
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Fig. 1. Images of comparable areas of breech face impressions on two NIST standard cartridge cases acquired by confocal 
microscopy after Gaussian filtering shown as (a) and (b). Fig. 1(c) is the transform of (b) using the algorithm proposed in 
Sec. 2 at the position of optimum registration, CCf,,,^^. = 97.9 %. 
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transformation consisting of x- and ^-translation and 
z-rotation. Furthermore, for an affine [2] transformation 
of images with different magnification scales,in case, 
for example, two images were measured with different 
lateral magnifications, this algorithm can be further 
expanded by adding a scaling factor /. The first and 
second order derivative matrix for implementing the 
optimization can then be derived from the following 
equation: 



x = T^ + l-{x- cos 9 - jsin 0) 
7 = 7, +l{xsmO + ycosO^ 



(15) 
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